/*

  xmunipack - fits image

  Copyright © 2009-2011 F.Hroch (hroch@physics.muni.cz)

  This file is part of Munipack.

  Munipack is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
  
  Munipack is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.

*/


#include "fits.h"
#include <cfloat>
#include <wx/wx.h>

#ifdef __WXDEBUG__
#include <wx/debug.h>
#include <wx/txtstrm.h>
#include <wx/wfstream.h>
#endif

using namespace std;

// ------------   FitsGeometry


FitsGeometry::FitsGeometry(const FitsArray& a): FitsArray(a) 
{ 
  wxASSERT(IsOk());
}

FitsArray FitsGeometry::Scale(float zoom)
{
  wxASSERT(zoom > 0.0);

  long *naxes = new long[Naxis()];
  for(int i = 0; i < Naxis(); i++) naxes[i] = max(int(zoom*Naxes(i)),1);
  FitsArray a = Scale(naxes);
  delete[] naxes;
  return FitsArray(a);
}

FitsArray FitsGeometry::Scale(const long *naxes)
{
  wxASSERT(naxes);
  long npixels = 1;
  for(int i = 0; i < Naxis(); i++) {
    wxASSERT(naxes[i] > 0);
    npixels = npixels*naxes[i];
  }

  float *a = new float[npixels];

  switch(Naxis()) {
    
  case 1:

    if( naxes[0] < Naxes(0) ) {
      // shrink

      float dx = max(double(Naxes(0))/double(naxes[0]),1.0);

      for(int i = 0; i < naxes[0]; i++) {
	float x = dx*(i + 0.5);
	*(a+i) = MeanLine(x,dx);
      }
    }
    else {
      // zoom

      float dx = float(Naxes(0))/float(naxes[0]);

      for(int i = 0; i < naxes[0]; i++) {
	int x = dx*(i + 0.5);
	*(a+i) = Pixel(x);
      }
    }
    break;

  case 2:

    if( naxes[0] < Naxes(0) && naxes[1] < Naxes(1) ) {
      // shrink

      float dx = max(double(Naxes(0))/double(naxes[0]),1.0);
      float dy = max(double(Naxes(1))/double(naxes[1]),1.0);

      for(int j = 0; j < naxes[1]; j++) {
	float y = dy*(j + 0.5);
	float *aa = a + j*naxes[0];

	for(int i = 0; i < naxes[0]; i++) {
	  float x = dx*(i + 0.5);
	  *(aa + i) = MeanRect(x,y,dx,dy);
	}
      }

    }
    else if( naxes[0] == Naxes(0) && naxes[1] == Naxes(1) ) {
      // 1:1

      FitsArrayData *data = static_cast<FitsArrayData *>(m_refData);
      wxASSERT(data);
      wxASSERT(data->array);
  
      int width = naxes[0];
      int height = naxes[1];
      int nbytes = width*sizeof(float);

      for(int i = 0; i < height; i++) {
	int n = i*width;
	memcpy(a+n,data->array+n,nbytes);
      }
    }
    else {
      // zoom

      float dx = float(Naxes(0))/float(naxes[0]);
      float dy = float(Naxes(1))/float(naxes[1]);

      for(int j = 0; j < naxes[1]; j++) {
	int y = dy*(j + 0.5);
	float *aa = a + j*naxes[0];

	for(int i = 0; i < naxes[0]; i++) {
	  int x = dx*(i + 0.5);
	  *(aa + i) = Pixel(x,y);
	}
      }
    }
    break;

  default:

    wxFAIL_MSG("FitsImageGray arbitrary dimension are not implemented yet");
    break;
  }

  long *ns = new long[Naxis()];
  for(int k = 0; k < Naxis(); k++)
    ns[k] = naxes[k];

  return FitsArray(*this,Naxis(),ns,a);
}


FitsArray FitsGeometry::Scale(int wsize, int hsize)
{
  wxASSERT(IsOk());
  wxASSERT(wsize > 0 || hsize > 0);

  if( Naxis() == 1 ) {

    const int step = 9;
    int w = Width();
    int h;
    
    int l = (wsize*hsize)/step;

    w = wsize;
    if( hsize > 0 ) 
      h = max(step*(hsize/step),1);
    else
      h = max(step*(wsize/step),1);
    l = w*(h/step);

    long ns[] = {Naxes(0)};
    FitsArray image = Scale(ns);

    float *a = new float[w*h];

    for(int i = 0; i < w*h; i++) a[i] = 0;

    for(int x = 0; x < l; x++) {
      int j = x/w;
      int i = x - j*w;
      float f = image.Pixel(x);
      for(int m = j*step;  m < (j + 1)*step-1; m++)
	a[m*w+i] = f;
    }
    long *nn = new long[2]; nn[0] = w; nn[1] = h;
    return FitsArray(*this,2,nn,a);
    
  }
  else if( Naxis() == 2 ) {

    float ratio;
    int w = Width();
    int h = Height();
  
    if( w > h )
      ratio = float(w)/float(wsize);
    else
      ratio = float(h)/float(hsize);

    w = max(int(w/ratio),1);
    h = max(int(h/ratio),1);

    long ns[] = {w,h};
    return Scale(ns);
  }

  wxFAIL_MSG("FitsGeom::Scale implemented for 1d and 2d images only");
  return FitsArray();
}

FitsArray FitsGeometry::SubArray(int xoff, int yoff, int w, int h)
{
  long npixels = w*h;
  long *ns = new long[2]; ns[0] = w; ns[1] = h;
  float *a = new float[npixels];
  int bytes;

  const FitsArrayData *data = static_cast<FitsArrayData *>(m_refData);
  wxASSERT(data);
  const float *array = data->array;
  wxASSERT(array);

  switch(Naxis()) {
    
  case 2:
    
    bytes = w*sizeof(float);
    for(int j = 0; j < h; j++)
      memcpy(a+j*w,array+(yoff+j)*data->naxes[0]+xoff,bytes);
    
    break;

  default:
    
    wxFAIL_MSG("FitsGeom::SubArray arbitrary dimension isn't implemented");
    break;
  }

  return FitsArray(*this,2,ns,a);  
}

void FitsGeometry::SetSubArray(int xoff, int yoff, const FitsArray& sub)
{
  // JUST ONLY *TWO* DIMENSIONS ARE IMPLEMENTED

  wxASSERT(sub.Naxis() == 2 && Naxis() == 2);
  FitsArrayData *data = static_cast<FitsArrayData *>(m_refData);
  wxASSERT(data);
  float *array = data->array;
  wxASSERT(array);

  wxASSERT(0 <= xoff && xoff + sub.Naxes(0) <= data->naxes[0]);
  wxASSERT(0 <= yoff && yoff + sub.Naxes(1) <= data->naxes[1]);

  int w = wxMin(sub.Naxes(0),data->naxes[0]-xoff);
  int h = wxMin(sub.Naxes(1),data->naxes[1]-yoff);

  /*
  if( w != sub.Naxes(0) || h != sub.Naxes(1) )
    wxLogDebug("%d %d %d %d %d %d %d %d", 
	       (int)h,(int)sub.Naxes(0),(int)w,(int) sub.Naxes(1),(int)data->naxes[0],(int)xoff,(int)data->naxes[1],(int)yoff);
  */

  const float *a = sub.PixelData();
  wxASSERT(a);
  int bytes = w*sizeof(float);

  for(int j = 0; j < h; j++)
    memcpy(array+(yoff+j)*data->naxes[0]+xoff,a+j*w,bytes);    
}

// ------------   FitsImage

FitsImage::FitsImage() {}

FitsImage::FitsImage(const vector<FitsArray>& a)
{
  wxASSERT(a.size() > 0);
  arrays = a;
}

FitsImage::FitsImage(const FitsArray& array)
{
  if( array.Naxis() == 3 )
    for(int l = 0; l < array.Naxes(2); l++)
      arrays.push_back(array.Plane(l));

  else
    arrays.push_back(array);
    

}

FitsImage::FitsImage(const FitsArray& r,const FitsArray& g,const FitsArray& b)
{
  arrays.push_back(r);
  arrays.push_back(g);
  arrays.push_back(b);
}

FitsImage::FitsImage(const FitsFile& fits, int sel)
{
  wxASSERT(fits.IsOk());

  if( sel == -1 && fits.Type() == FITS_COLOR ) {
    // color 

    for(int k = 1; k <= 3; k++) {
      FitsHdu hdu = fits.Hdu(k);
      arrays.push_back(FitsArray(hdu));
    }
  }
  else if( 0<=sel&&sel<int(fits.HduCount()) && fits.Hdu(sel).Type()==HDU_IMAGE){
    // gray
 
    wxASSERT( fits.Hdu(sel).Type() == HDU_IMAGE );
    arrays.push_back(fits.Hdu(sel));
  }    
  else
    wxFAIL_MSG("FitsImage::FitsImage"+fits.GetName());
}

FitsImage::~FitsImage() 
{ 
  arrays.clear();
}

bool FitsImage::IsOk() const 
{ 
  return arrays.size() > 0;
}

bool FitsImage::IsColor() const 
{ 
  return arrays.size() > 1;
}

const std::vector<FitsArray> FitsImage::GetArrays() const
{
  return arrays;
}


int  FitsImage::GetWidth() const 
{
  if( arrays.empty() )
    return 0;
  else
    return arrays[0].Width();
}

int  FitsImage::GetHeight() const 
{
  if( arrays.empty() )
    return 0;
  else
    return arrays[0].Height();
}

size_t FitsImage::GetCount() const
{
  return arrays.size();
}

FitsArray FitsImage::Item(size_t k) const
{ 
  wxASSERT(0 <= k && k < arrays.size());
  return arrays[k];
}

FitsImage FitsImage::Scale(float zoom)
{    
  vector<FitsArray> temp;

  for(vector<FitsArray>::const_iterator k=arrays.begin(); k!=arrays.end();++k){
    FitsGeometry g(*k);
    temp.push_back(g.Scale(zoom));
  }

  return FitsImage(temp);
}

FitsImage FitsImage::Scale(int width, int height)
{    
  vector<FitsArray> temp;

  for(vector<FitsArray>::const_iterator k=arrays.begin(); k!=arrays.end();++k){
    FitsGeometry g(*k);
    temp.push_back(g.Scale(width,height));
  }

  return FitsImage(temp);
}

void FitsImage::Rescale(float zoom)
{
  vector<FitsArray> temp;

  wxASSERT(zoom > 0.0);

  for(vector<FitsArray>::const_iterator k=arrays.begin(); k!=arrays.end();++k){
    FitsHdu hdu(*k);
    int naxis = k->Naxis();
    long *naxes = new long[naxis];
    for(int i = 0; i < naxis; i++) naxes[i] = max(int(zoom*k->Naxes(i)),1);
    long n = 1;
    for(int i = 0; i < naxis; i++) n = n*naxes[i];
    float *a = new float[n];
    temp.push_back(FitsArray(hdu,naxis,naxes,a));
  }

  arrays = temp;
}

FitsImage FitsImage::GetSubImage(int x0, int y0, int w, int h)
{
  vector<FitsArray> temp;

  for(vector<FitsArray>::const_iterator k=arrays.begin(); k!=arrays.end();++k){
    wxASSERT(k->Naxis() == 2);
    FitsGeometry g(*k);
    temp.push_back(g.SubArray(x0,y0,w,h));
  }

  return FitsImage(temp);


  // if( arrays.size() == 1 ) {

  //   FitsGeometry g(arrays[0]);

  //   if( g.Naxis() != 2 )
  //     wxFAIL_MSG("FitsImage::SubImage() only 2D images implemented yet");

  //   return FitsImage(g.SubArray(x0,y0,w,h));
  // }
  // else if( arrays.size() == 3 ) {

  //   FitsGeometry r(arrays[0]);
  //   FitsGeometry g(arrays[1]);
  //   FitsGeometry b(arrays[2]);

  //   return FitsImage(r.SubArray(x0,y0,w,h),g.SubArray(x0,y0,w,h),
  // 		     b.SubArray(x0,y0,w,h));
  // }

  // wxFAIL_MSG("FitsImage::SubImage only gray or color images implemented");
  // return FitsImage();
}


FitsImage FitsImage::Thumb(int wsize, int hsize)
{
  return Scale(wsize,hsize);
}      


void FitsImage::SetSubImage(int x, int y, const FitsImage& srcimage)
{
  vector<FitsArray> src(srcimage.GetArrays());

  wxASSERT(srcimage.GetCount() == GetCount());

  for(size_t k = 0; k < arrays.size(); k++) {
    FitsGeometry g(arrays[k]);
    g.SetSubArray(x,y,src[k]);
  }
}


// ---- FitsDisplay  ------

FitsDisplay::FitsDisplay(const wxString& t): 
  cspace(t),temp(6500.0),xtemp(1.0/0.950456),ytemp(1.0),ztemp(1.0/1.088754) 
{
  spaces.Add("sRGB");
  spaces.Add("AdobeRGB");

  if( t.IsEmpty() )
    cspace = spaces[0];

  wxASSERT(spaces.Index(cspace) != wxNOT_FOUND);
}

FitsPalette FitsDisplay::GetPalette() const
{
  return pal;
}

FitsItt FitsDisplay::GetItt() const
{
  return itt;
}

FitsColor FitsDisplay::GetColor() const
{
  return color;
}

float FitsDisplay::GetTemperature() const
{
  return temp;
}

void FitsDisplay::SetItt(const FitsItt& i)
{
  itt = i;
}

void FitsDisplay::SetPalette(const FitsPalette& p)
{
  pal = p;
}

void FitsDisplay::SetColor(const FitsColor& c)
{
  color = c;
}

void FitsDisplay::SetTemperature(float t)
{
  temp = t;

  // http://en.wikipedia.org/wiki/Planckian_locus
  // we compute factors to for color balance of XYZ values
  // all computations are available only for 4000 < T < 25000 K
  // the used approximation looks unaccurate

  double t2 = t*t;
  double t3 = t*t2;

  double xc = 0.24039 + 0.2226347e3/t + 2.1070379e6/t2 - 3.0258469e9/t3;
  double x2 = xc*xc;
  double x3 = xc*x2;
  double yc = 3.0817580*x3 - 5.87338670*x2 + 3.75112997*xc - 0.37001483;
  double zc = 1.0 - xc - yc;

  double tx = xc/yc;
  double tz = zc/yc;

  xtemp = 1.0/tx;
  ztemp = 1.0/tz;

  //  wxLogDebug(_("%f %f %f %f %f"),xtemp,ztemp,xc,yc,zc);
}

void FitsDisplay::SetColorspace(const wxString& cs)
{
  wxASSERT(spaces.Index(cs) != wxNOT_FOUND);
  cspace = cs;
}

wxString FitsDisplay::GetColorspace() const
{
  return cspace;
}


wxArrayString FitsDisplay::GetArraySpaces() const
{
  return spaces;
}


FitsBitmap FitsDisplay::ConvertOfGray(const FitsArray& array)
{
  // This part does gray scale image transformation to 
  // to an output space.
  //
  // For gray images, only gamma-like-correction is required.

  wxASSERT(array.IsOk());

  int npix = array.Npixels();

  float *f = itt.Scale(npix,array.PixelData());

  float *g = 0;
  if( cspace == "sRGB" )
    g = sRGBGamma(npix,f);
  else if( cspace == "AdobeRGB" )
    g = AdobeGamma(npix,f);
  else
    wxFAIL_MSG("An unknown display colorspace");
  delete[] f;

  unsigned char *rgb = pal.RGB(npix,g);

  delete[] g;

  return FitsBitmap(array.Width(),array.Height(),rgb);
}

float *FitsDisplay::sRGBGamma(int n, float *a)
{
  float *f = new float[n];
  for(int i = 0; i < n; i++) {
    float r = a[i];
    if( r < 0.00304f )
      f[i] = 12.92f*r;
    else
      f[i] = 1.055f*powf(r,0.4166667f) - 0.055f;
  }
  return f;
}


float *FitsDisplay::AdobeGamma(int n, float *a)
{
  const float q = 1.0/2.19921875;

  float *f = new float[n];
  for(int i = 0; i < n; i++)
    f[i] = powf(a[i],q);
  return f;
}


// FitsBitmap FitsDisplay::ConvertOfGray(const FitsArray& array)
// {
//   // This part does gray scale image transformation to 
//   // to an output space.
//   //
//   // For gray images, only gamma-correction is required.

//   wxASSERT(array.IsOk());

//   int width = array.Width();
//   int height = array.Height();
//   int npix = width*height;

//   unsigned char *rgb = (unsigned char *) malloc(3*npix);
//   unsigned char *p = rgb;
//   const float *a = array.PixelData();
//   const unsigned char *rpal = pal.RData();
//   const unsigned char *gpal = pal.GData();
//   const unsigned char *bpal = pal.BData();
//   int npal = pal.GetColors() - 1;

//   for(int i = 0; i < npix; i++) {
//     float x = Gamma(itt.Scale(a[i]));
//     //float x = color.Gamma(itt.Scale(array.Pixel(i)));

//     int l = int(255.0*(x + 1.0));
//     if( l < 0 )
//       l = 0;
//     else if( l > npal )
//       l = npal;

//     *p++ = rpal[l];
//     *p++ = gpal[l];
//     *p++ = bpal[l];

//     /*
//     unsigned char r,g,b;
//     pal.RGB(x,r,g,b);
//     *p++ = r;
//     *p++ = g;
//     *p++ = b;
//     */
//   }
//   return FitsBitmap(width, height, rgb);
// }

// FitsBitmap FitsDisplay::ConvertOfColor(const FitsArray& red,
// 				       const FitsArray& green,
// 				       const FitsArray& blue)
// {
//   // This part does color transformation from an input color space 
//   // to an output space.
  
//   // Digital cameras usually produces colors in sRGB or AdobeRGB.
  
//   // Astronomical filter series usually stores data in their proper systems.
//   // Plain usage leads to a deformation of colors.
//   // Documents 
//   //         http://www.fho-emden.de/~hoffmann/ciexyz29082000.pdf
//   //     and http://www.fho-emden.de/~hoffmann/cielab03022003.pdf 
//  //          http://www.babelcolor.com/download/A%20review%20of%20RGB%20color%20spaces.pdf
//   // detaily describes the general transformation.

//   // Therefore, we are suppose that the input is in XYZ color space.
//   // (XYZ should be obtained from any RAW digital photography with dcraw).
//   // (XYZ should be obtained from Landolt's BVR by a custom transformation.)

//   // http://en.wikipedia.org/wiki/RGB_color_spaces
//   // http://en.wikipedia.org/wiki/CIE_1931_color_space
//   // http://en.wikipedia.org/wiki/D65
  
//   // The main purpose of this routine is scale luminosity with
//   // respecting of colors. This is done with help of xyY color system. 

//   // Luminosity scaled XYZ is directly transformed to sRGB.
//   // (Adobe RGB (Mac) or sRGB is not directly supported by wxWidgets.)

//   wxASSERT(red.IsOk() && green.IsOk() && blue.IsOk());
//   wxASSERT(red.Width() == green.Width() && red.Width() == blue.Width());
//   wxASSERT(red.Height() == green.Height() && red.Height() == blue.Height());

//   int width = red.Width();
//   int height = red.Height();
//   int npix = width*height;

//   int ncolors = color.GetColors();
//   int nbands = color.GetBands();
//   //  float *cb = new float[ncolors*nbands];
//   float cb[ncolors][nbands];
//   for(int i = 0; i < ncolors; i++)
//     for(int j = 0; j < nbands; j++)
//       cb[i][j] = color.GetTrans(i,j);

//   float weight[nbands], level[nbands];
//   for(int j = 0; j < nbands; j++) {
//     weight[j] = color.GetWeight(j);
//     level[j] = color.GetLevel(j);
//   }

//   // The observed values are supposed as XYZ (not xyY) without scaling
//   // into interval 0 to 100. There, we set its white point to D65 and
//   // color weights by user or camera:

//   // There, we are converting the wider dynamic range to a small one.
//   // To preserve colors and speed up computation, we scale only
//   // luminosity (by defined human's eye) and use 'color indexes'.

//   unsigned char *rgb = (unsigned char *) malloc(3*npix);
//   unsigned char *p = rgb;
//   float f = itt.InvScale(1/*255*/);
//   float z0 = itt.InvScale(0);
//   f = 1.0/itt.GetSensitivity();
//   z0 = itt.GetBlack();
//   float hue = color.GetHue();
//   float satur = color.GetSaturation();
//   bool nvision = color.GetNightVision();
//   //  float t = color.GetNightThresh();
//   //  float w = color.GetNightWidth();

//   //  wxLogDebug(wxT("%f %f %f %f"),f,z,itt.GetMad(),itt.GetMed());

//   //  wxLogDebug(wxT("%f %f"),color.GetWhitePoint_U(),color.GetWhitePoint_V());

//   FitsItt litt(itt.GetItt());
//   //  litt.SetOffset(50.0);
//   //  litt.SetOffset(0.0);
//   //litt.SetSlope(0.05/0.04);

//   /*
//   const float *x = red.PixelData();
//   const float *y = green.PixelData();
//   const float *z = blue.PixelData();
//   */
    
//   const float *x = blue.PixelData();
//   const float *y = green.PixelData();
//   const float *z = red.PixelData();
  

//   for(int i = 0; i < npix; i++) {

//     float X = x[i];
//     float Y = y[i];
//     float Z = z[i];

//     /*    
//     float X = red.Pixel(i);
//     float Y = green.Pixel(i);
//     float Z = blue.Pixel(i);
//     */

//     //    float grass = green.Pixel(i);
//     //    float f = itt.Scale(Y);

//     //wxLogDebug(wxT("%f %f"),f,itt.InvScale(255));

//     // implementovat roztahovani slideru na intenzitu podle dat z histogramu
//     // right click v obrazku aktivuje menu s copy na kopirovani pixelu

//     //    float z = itt.InvScale(0);


//     //    float B = X;
//     //    float V = Y;
//     //    float R = Z;

//     // rotate by color matrix
//     //    float d[3] = {Z-20700.0,Y-9500.0,(X-3780.0)*1.0};
//     //    float d[3] = {(Z-3830.0)*0.6,Y-8780.0,X-19700.0};
//     //float d[3] = {Z,Y,X};
//     //    float d[3] = {0.6*(B-3778.0),V-9478.0,R-20740.0};
//     //    float d[3] = {0.6*(B-3855.0),V-8785.0,R-19210.0};
//     //    float d[3] = {B,V,R};
//     float d[3] = {X,Y,Z};
//     /*
//     for(int i = 0; i < nbands; i++)
//       d[i] = weight[i]*(d[i] - level[i]);
//     */
//       //      d[i] = color.GetWeight(i)*(d[i] - color.GetLevel(i));

//     //    float t[ncolors];
    
//     float *c[] = { &Z, &Y, &X };
//     for(int i = 0; i < ncolors; i++) {
//       float s = 0.0;
//       for(int j = 0; j < nbands; j++)
// 	s = s + cb[j][i]*d[j];
//       //t[i] = s;
//       **(c+i) = s;
//     }
    
//     /*
//     X = t[2];
//     Y = t[1];
//     Z = t[0];
//     */

//     //    color.Transform(3,d,X,Y,Z);

//     //    if( i == 2000 )
//     //      wxLogDebug(_("%f %f %f"),X,Y,Z);

//     // ???
    
//     X = X > z0 ? X - z0 : 0.0;
//     Y = Y > z0 ? Y - z0 : 0.0;
//     Z = Z > z0 ? Z - z0 : 0.0;
    

//     /*
//     if( X > 100.0 ) X = 100.0;
//     if( Y > 100.0 ) Y = 100.0;
//     if( Z > 100.0 ) Z = 100.0;
//     */

//     /*
//     X = X - z;
//     Y = Y - z;
//     Z = Z - z;
//     */

//     float L,u,v;
//     color.XYZ_Luv(X,Y,Z,f,L,u,v);

//     //    float LL = L;

//     L = 100.0*litt.Scale(L/100.0);
//     //    L = 100.0*litt.Scale((L-50.0)*0.05);
//     //L = 100.0/(1.0+exp(-(L-50)*0.1));
//     //    L = 100.0*(erf((L-50.0)*0.05)+1.0);
//     //itt.Scale(f*L+z);//itt.InvScale(2.55*L)/f*100.0;
//     //    color.XYZ_Luv(X-z,Y-z,Z-z,/*Y/f*/f-z,L,u,v);
//     /*
//     if( L > 100.0 ) u = v = 0.0;
//     if( L > 100.0 ) L = 100.0;
//     if( L < 0.0 ) L = 0.0;
//     */

//     // color saturation
//     /*
//     float hue = atan2f(v,u);
//     float chroma = sqrtf(v*v + u*u);
//     u = color.GetSaturation()*chroma*cosf(hue);
//     v = color.GetSaturation()*chroma*sinf(hue);

//     if( color.GetNightVision() ) {
//       float t = color.GetNightThresh();
//       float w = color.GetNightWidth();
 
//       float p = 1.0/(1.0 + exp(-(Y - z - t)/w));
//       u = p*u;
//       v = p*v;
//     }
//     color.Luv_XYZ(L,u,v,1.0,X,Y,Z);
//     */

//     // color adjustment
//     float h = (v != 0.0 && u != 0.0 ? atan2f(v,u) : 0.0) + hue;
//     float ch = satur*sqrtf(v*v + u*u);
//     u = ch*cosf(h);
//     v = ch*sinf(h);
    
//     color.Luv_XYZ(L,u,v,1.0,X,Y,Z);
    
//     if( nvision ) {
      
//       //      color.Luv_XYZ(L,u,v,1.0,X,Y,Z);
//       float w = color.NightProfile(L);
//       float s = color.Scotopic(X,Y,Z);

//       float X1,Y1,Z1;
//       //      color.Luv_XYZ(100.0*s,0.0,0.0,1.0,X1,Y1,Z1);
//       // simplified transformation for gray, 
//       //   just only temperature correction is needed
//       X1 = s;//*0.950456;
//       Y1 = s;
//       Z1 = s;//*1.088754;

//       //      color.XYZ_Luv(X,Y,Z,f,L,u,v);

//       //u = w*u;
//       //      v = w*v;
//       //      L = s;//w*L + (1.0 - w)*s;

//       float w1 = 1.0 - w;
//       X = w*X + w1*X1;
//       Y = w*Y + w1*Y1;
//       Z = w*Z + w1*Z1;
//     }
    
//     //    color.Luv_XYZ(L,u,v,1.0,X,Y,Z);
    
//     //    OutputRGB(X,Y,Z,R,G,B);

//     unsigned char r,g,b;
//     color.XYZ_sRGB(X,Y,Z,r,g,b);
//     //    color.XYZ_sRGB(/*flux,*/red.Pixel(i),grass,blue.Pixel(i),r,g,b);
//     *p++ = r;
//     *p++ = g;
//     *p++ = b;

//     /*
//     color.ConvertRGB(X,Y,Z,r,g,b);
//     *p++ = r;
//     *p++ = g;
//     *p++ = b;    
//     */
//   }
//   //  delete[] cb;
//   return FitsBitmap(width, height, rgb);
// }




FitsBitmap FitsDisplay::ConvertOfColorX(const std::vector<FitsArray>& band)
{
  // This part does color transformation from an instrumental color space 
  // to a display color space. Optionaly tune color parameters.
  
  // Digital cameras usually produces colors in sRGB or AdobeRGB.
  
  // Astronomical filter series usually stores data in their proper systems.
  // Plain usage leads to a deformation of colors.
  // Documents 
  //         http://www.fho-emden.de/~hoffmann/ciexyz29082000.pdf
  //     and http://www.fho-emden.de/~hoffmann/cielab03022003.pdf 
  // detaily describes the general transformation.

  // Therefore, we are converting the input to the XYZ color space.
  // (XYZ should be obtained from any RAW digital photography with dcraw).
  // (XYZ should be obtained from Landolt's BVR by a custom transformation.)

  // http://en.wikipedia.org/wiki/RGB_color_spaces
  // http://en.wikipedia.org/wiki/CIE_1931_color_space
  // http://en.wikipedia.org/wiki/D65
  
  // Luminosity scaled XYZ is transformed to some RGB.
  // (Adobe RGB (Mac) or sRGB is not directly supported by wxWidgets.)

#ifdef __WXDEBUG__
  wxASSERT(band.size() > 0);
  wxASSERT(band.size() == (size_t) color.GetBands());
  for(size_t i = 0; i < band.size(); i++) {
    wxASSERT(band[i].IsOk());
    if( i > 0 ) {
      wxASSERT(band[i].Naxis() == band[i-1].Naxis());
      for(size_t j = 0; j < (size_t) band[i].Naxis(); j++)
	wxASSERT(band[i].Naxes(j) == band[i-1].Naxes(j));
    }
  }
#endif  

  wxASSERT(band[0].Npixels() == band[1].Npixels());
  wxASSERT(band[1].Npixels() == band[2].Npixels());

  int width = band[0].Width();
  int height = band[0].Height();
  int npix = band[0].Npixels();
  float *X = new float[npix];
  float *Y = new float[npix];
  float *Z = new float[npix];

  // pointers to raw data
  const float **d = new const float*[band.size()];
  for(size_t i = 0; i < band.size(); i++)
    *(d+i) = band[i].PixelData();


  //  wxLogDebug("%d %d",(int)band.size(),(int)npix);


  // Instrumental to XYZ
  color.Instr_XYZ(npix,band.size(),d,X,Y,Z);

  delete[] d;
  /*
  for(size_t i = 0; i < band.size(); i++)
    delete[] *(d+i);
  */

  float *L = new float[npix];
  float *u = new float[npix];
  float *v = new float[npix];
  
  FitsItt litt(itt.GetItt());
  litt.Init(0.0,1.0);
  litt.SetKind(ITT_KIND_ABS);
  litt.SetBlack(0.0);
  litt.SetContrast(1.0);
  litt.SetAmp(itt.GetAmp());
  litt.SetZero(itt.GetZero());
  
  float f = itt.InvScale(1);
  float z0 = itt.InvScale(0);
  //  f = 1.0/itt.GetSensitivity();
  f = itt.GetContrast();
  z0 = itt.GetBlack();
  //wxLogDebug("%f %f",f,z0);
  //  litt.SetBlack(z0);
  
  // prescale
  for(int i = 0; i < npix; i++) {

    float t = X[i] - z0;
    X[i] = t > 0 ? t : 0.0;

    t = Y[i] - z0;
    Y[i] = t > 0.0 ? t : 0.0;

    t = Z[i] - z0;
    Z[i] = t > 0.0 ? t : 0.0;
  }
  
  // XYZ to Luv
  color.XYZ_Luv(npix,X,Y,Z,f,L,u,v);

  // scale luminosity
  for(int i = 0; i < npix; i++) L[i] = L[i]/100.0;
  float *xL = litt.Scale(npix,L);
  for(int i = 0; i < npix; i++) L[i] = xL[i]*100.0;
  delete[] xL;

  // tune saturation and hue
  color.TuneColors(npix,L,u,v);

  // convert back to XYZ
  color.Luv_XYZ(npix,L,u,v,1.0,X,Y,Z);

  // night vision
  color.NightVision(npix,L,X,Y,Z);

  // disply colorspace
  unsigned char *rgb = 0;

  if( cspace == "sRGB" )
    rgb = XYZ_sRGB(npix,X,Y,Z);
  else if( cspace == "AdobeRGB" )
    rgb = XYZ_AdobeRGB(npix,X,Y,Z);
  else
    wxFAIL_MSG("An unknown display colorspace");

  wxASSERT(rgb);

  delete[] X;
  delete[] Y;
  delete[] Z;
  delete[] L;
  delete[] u;
  delete[] v;

  return FitsBitmap(width,height,rgb);  
}



//  Generalization for multiple colors
FitsBitmap FitsDisplay::ConvertOfColor(const std::vector<FitsArray>& band)
{
  // This part does color transformation from an input color space 
  // to an output space.
  
  // Digital cameras usually produces colors in sRGB or AdobeRGB.
  
  // Astronomical filter series usually stores data in their proper systems.
  // Plain usage leads to a deformation of colors.
  //
  // Documents 
  //         http://www.fho-emden.de/~hoffmann/ciexyz29082000.pdf
  //     and http://www.fho-emden.de/~hoffmann/cielab03022003.pdf 
  // detaily describes the general transformation.

  // Therefore, we are suppose that the input is in XYZ color space.
  // (XYZ should be obtained from any RAW digital photography with dcraw).
  // (XYZ should be obtained from Landolt's BVR by a custom transformation.)

  // http://en.wikipedia.org/wiki/RGB_color_spaces
  // http://en.wikipedia.org/wiki/CIE_1931_color_space
  // http://en.wikipedia.org/wiki/D65
  
  // The main purpose of this routine is scale luminosity with
  // respecting of colors. This is done with help of xyY color system. 

  // Luminosity scaled XYZ is directly transformed to sRGB.
  // (Adobe RGB (Mac) or sRGB is not directly supported by wxWidgets.)

#ifdef __WXDEBUG__
  wxASSERT(band.size() > 0);
  wxASSERT(band.size() == (size_t) color.GetBands());
  for(size_t i = 0; i < band.size(); i++) {
    wxASSERT(band[i].IsOk());
    if( i > 0 ) {
      wxASSERT(band[i].Naxis() == band[i-1].Naxis());
      for(size_t j = 0; j < (size_t) band[i].Naxis(); j++)
	wxASSERT(band[i].Naxes(j) == band[i-1].Naxes(j));
    }
  }
#endif  

  int width = band[0].Width();
  int height = band[0].Height();
  //  int npix = width*height;
  int npix = band[0].Npixels();

  int ncolors = color.GetColors();
  int nbands = color.GetBands();
  //  float *cb = new float[ncolors*nbands];
  float cb[ncolors][nbands];
  for(int i = 0; i < ncolors; i++)
    for(int j = 0; j < nbands; j++) {
      cb[i][j] = color.GetTrans(i,j);
      //      wxLogDebug(_("%d %d %f"),i,j,cb[i][j]);
    }

  float weight[nbands], level[nbands];
  for(int j = 0; j < nbands; j++) {
    weight[j] = color.GetWeight(j);
    level[j] = color.GetLevel(j);
  }

  // The observed values are supposed as XYZ (not xyY) without scaling
  // into interval 0 to 100. There, we set its white point to D65 and
  // color weights by user or camera:

  // There, we are converting the wider dynamic range to a small one.
  // To preserve colors and speed up computation, we scale only
  // luminosity (by defined human's eye) and use 'color indexes'.

  //  unsigned char *rgb = (unsigned char *) malloc(3*npix);
  //  unsigned char *p = rgb;
  float f = itt.InvScale(1/*255*/);
  float z0 = itt.InvScale(0);
  f = 1.0/itt.GetSensitivity();
  z0 = itt.GetBlack();
  float hue = color.GetHue();
  float satur = color.GetSaturation();
  bool nvision = color.GetNightVision();
  //  float t = color.GetNightThresh();
  //  float w = color.GetNightWidth();

  //  wxLogDebug(wxT("%f %f %f %f"),f,z,itt.GetMad(),itt.GetMed());

  //  wxLogDebug(wxT("%f %f"),color.GetWhitePoint_U(),color.GetWhitePoint_V());

  FitsItt litt(itt.GetItt());
  //  litt.SetOffset(50.0);
  //  litt.SetOffset(0.0);
  //litt.SetSlope(0.05/0.04);

  /*
  const float *x = red.PixelData();
  const float *y = green.PixelData();
  const float *z = blue.PixelData();
  */

  /*    
  const float *x = blue.PixelData();
  const float *y = green.PixelData();
  const float *z = red.PixelData();
  */
  //  float **d = (float **) malloc(band.size()*sizeof(float *));//new *float[band.size()];
  const float **d = new const float*[band.size()];
  for(size_t i = 0; i < band.size(); i++)
    *(d+i) = band[i].PixelData();


  float X,Y,Z;
  float *c[] = { &X, &Y, &Z };
  //  float *c[] = { &Z, &Y, &X };


  float x[npix],y[npix],z[npix];


  if( color.GetColorspace().Find("XYZ") != wxNOT_FOUND ) {
    wxASSERT(ncolors == nbands);

    memcpy(x,band[2].PixelData(),npix*sizeof(float));
    memcpy(y,band[1].PixelData(),npix*sizeof(float));
    memcpy(z,band[0].PixelData(),npix*sizeof(float));
  }
  else {

    float b[nbands];
    for(int n = 0; n < npix; n++) {

      for(int j = 0; j < nbands; j++) {
	b[j] = weight[j]*(d[j][n] - level[j]);
	if( b[j] < 0.0 ) b[j] = 0.0;
      }

      for(int i = 0; i < ncolors; i++) {
	float s = 0.0;
	for(int j = 0; j < nbands; j++)
	  s = s + cb[i][j]*b[j];
	**(c+i) = s;
	//	**(c+(ncolors-i-1)) = s;
      }

      x[n] = *c[2];
      y[n] = *c[1];
      z[n] = *c[0];

      /*
      x[n] = b[0];
      y[n] = b[1];
      z[n] = b[2];
      */
    }
  }

  //  color.XYZ_Luv(npix,x,y,z,f,L,u,v);



  for(int n = 0; n < npix; n++) {

    /*
    float X = x[i];
    float Y = y[i];
    float Z = z[i];
    */

    /*    
    float X = red.Pixel(i);
    float Y = green.Pixel(i);
    float Z = blue.Pixel(i);
    */

    //    float grass = green.Pixel(i);
    //    float f = itt.Scale(Y);

    //wxLogDebug(wxT("%f %f"),f,itt.InvScale(255));

    // implementovat roztahovani slideru na intenzitu podle dat z histogramu
    // right click v obrazku aktivuje menu s copy na kopirovani pixelu

    //    float z = itt.InvScale(0);


    //    float B = X;
    //    float V = Y;
    //    float R = Z;

    // rotate by color matrix
    //    float d[3] = {Z-20700.0,Y-9500.0,(X-3780.0)*1.0};
    //    float d[3] = {(Z-3830.0)*0.6,Y-8780.0,X-19700.0};
    //float d[3] = {Z,Y,X};
    //    float d[3] = {0.6*(B-3778.0),V-9478.0,R-20740.0};
    //    float d[3] = {0.6*(B-3855.0),V-8785.0,R-19210.0};
    //    float d[3] = {B,V,R};
    //    float d[3] = {X,Y,Z};
    /*
    for(int i = 0; i < nbands; i++)
      d[i][n] = weight[i]*(d[i][n] - level[i]);
    */
      //      d[i] = color.GetWeight(i)*(d[i] - color.GetLevel(i));

    //    float t[ncolors];
    
    //    float *c[] = { &Z, &Y, &X };

//     wxString t = color.GetColorspace();
//     if( t.IsSameAs(wxT("XYZ")) ) {

//       wxASSERT(ncolors == nbands);
//       for(int i = 0; i < ncolors; i++)
// 	//**(c+i) = d[i][n];
// 	**(c+i) = d[ncolors-i-1][n];

//     } else {

//     for(int i = 0; i < ncolors; i++) {
//       float s = 0.0;
//       for(int j = 0; j < nbands; j++) {
// 	float w = weight[j]*(d[j][n] -/*band[j].Pixel(n)*/level[j]);
// 	s = s + cb[i][j]* (w > 0.0 ? w : 0.0);
//       }
//       //t[i] = s;
//       **(c+(ncolors-i-1)) = s;
//       //      **(c+i) = s;
//     }
//     }
    
    /*
    X = t[2];
    Y = t[1];
    Z = t[0];
    */

    //    color.Transform(3,d,X,Y,Z);

    //    if( i == 2000 )
    //      wxLogDebug(_("%f %f %f"),X,Y,Z);

    // ???

    X = x[n];
    Y = y[n];
    Z = z[n];
    
    X = X > z0 ? X - z0 : 0.0;
    Y = Y > z0 ? Y - z0 : 0.0;
    Z = Z > z0 ? Z - z0 : 0.0;
    

    /*
    if( X > 100.0 ) X = 100.0;
    if( Y > 100.0 ) Y = 100.0;
    if( Z > 100.0 ) Z = 100.0;
    */

    /*
    X = X - z;
    Y = Y - z;
    Z = Z - z;
    */

    float L,u,v;
    color.XYZ_Luv(X,Y,Z,f,L,u,v);

    //    float LL = L;

    L = 100.0*litt.Scale(L/100.0);

    //    L = 100.0*litt.Scale((L-50.0)*0.05);
    //L = 100.0/(1.0+exp(-(L-50)*0.1));
    //    L = 100.0*(erf((L-50.0)*0.05)+1.0);
    //itt.Scale(f*L+z);//itt.InvScale(2.55*L)/f*100.0;
    //    color.XYZ_Luv(X-z,Y-z,Z-z,/*Y/f*/f-z,L,u,v);
    /*
    if( L > 100.0 ) u = v = 0.0;
    if( L > 100.0 ) L = 100.0;
    if( L < 0.0 ) L = 0.0;
    */

    // color saturation
    /*
    float hue = atan2f(v,u);
    float chroma = sqrtf(v*v + u*u);
    u = color.GetSaturation()*chroma*cosf(hue);
    v = color.GetSaturation()*chroma*sinf(hue);

    if( color.GetNightVision() ) {
      float t = color.GetNightThresh();
      float w = color.GetNightWidth();
 
      float p = 1.0/(1.0 + exp(-(Y - z - t)/w));
      u = p*u;
      v = p*v;
    }
    color.Luv_XYZ(L,u,v,1.0,X,Y,Z);
    */

    // color adjustment
    float h = (v != 0.0 && u != 0.0 ? atan2f(v,u) : 0.0) + hue;
    float ch = satur*sqrtf(v*v + u*u);
    u = ch*cosf(h);
    v = ch*sinf(h);
    
    color.Luv_XYZ(L,u,v,1.0,X,Y,Z);
    
    if( nvision ) {
      
      //      color.Luv_XYZ(L,u,v,1.0,X,Y,Z);
      float w = color.NightProfile(L);
      float s = color.Scotopic(X,Y,Z);

      float X1,Y1,Z1;
      //      color.Luv_XYZ(100.0*s,0.0,0.0,1.0,X1,Y1,Z1);
      // simplified transformation for gray, 
      //   just only temperature correction is needed
      X1 = s;//*0.950456;
      Y1 = s;
      Z1 = s;//*1.088754;

      //      color.XYZ_Luv(X,Y,Z,f,L,u,v);

      //u = w*u;
      //      v = w*v;
      //      L = s;//w*L + (1.0 - w)*s;

      float w1 = 1.0 - w;
      X = w*X + w1*X1;
      Y = w*Y + w1*Y1;
      Z = w*Z + w1*Z1;
    }
    
    //    color.Luv_XYZ(L,u,v,1.0,X,Y,Z);
    
    //    OutputRGB(X,Y,Z,R,G,B);

    x[n] = X;
    y[n] = Y;
    z[n] = Z;

//     unsigned char r,g,b;
//     color.XYZ_sRGB(X,Y,Z,r,g,b);
//     //    color.XYZ_sRGB(/*flux,*/red.Pixel(i),grass,blue.Pixel(i),r,g,b);
//     *p++ = r;
//     *p++ = g;
//     *p++ = b;

    /*
    color.ConvertRGB(X,Y,Z,r,g,b);
    *p++ = r;
    *p++ = g;
    *p++ = b;    
    */
  }
  //  delete[] cb;


  unsigned char *rgb = 0;

  if( cspace == "sRGB" )
    rgb = XYZ_sRGB(npix,x,y,z);
  else if( cspace == "AdobeRGB" )
    rgb = XYZ_AdobeRGB(npix,x,y,z);
  else
    wxFAIL_MSG("An unknown display colorspace");

  wxASSERT(rgb);

  return FitsBitmap(width,height,rgb);
}


FitsBitmap FitsDisplay::GetImage(const FitsImage& image)
{
  if( image.IsColor() ) 
    return ConvertOfColorX(image.GetArrays());
  else
    return ConvertOfGray(image.Item(0));
}

unsigned char *FitsDisplay::Crange(size_t n, const float *a)
{
  unsigned char *f = new unsigned char[n];

  for(size_t i = 0; i < n; i++) {

    int j = int(255.0f * *a++);

    if( j < 0 )
      j = 0;
    else if( j > 255 )
      j = 255;

    f[i] = j;
  }

  return f;
}

unsigned char *FitsDisplay::XYZ_sRGB(size_t n, const float *X, const float *Y,
				     const float *Z)
{
  // http://en.wikipedia.org/wiki/SRGB

  float *r = new float[n];
  float *g = new float[n];
  float *b = new float[n];

  for(size_t i = 0; i < n; i++) {
    float x = X[i]/xtemp;
    float y = Y[i]/ytemp;
    float z = Z[i]/ztemp;

    r[i] =  3.2410f*x - 1.5374f*y - 0.4986f*z;
    g[i] = -0.9692f*x + 1.8760f*y + 0.0416f*z;
    b[i] =  0.0556f*x - 0.2040f*y + 1.0570f*z;
  }

  wxASSERT(cspace == "sRGB");

  float *r1 = sRGBGamma(n,r);
  float *g1 = sRGBGamma(n,g);
  float *b1 = sRGBGamma(n,b);

  delete[] r;
  delete[] g;
  delete[] b;

  unsigned char *ur = 0;
  unsigned char *ug = 0;
  unsigned char *ub = 0;

  ur = Crange(n,r1);
  ug = Crange(n,g1);
  ub = Crange(n,b1);

  delete[] r1;
  delete[] g1;
  delete[] b1;


  unsigned char *rgb = (unsigned char *) malloc(3*n);
  unsigned char *p = rgb;

  for(size_t i = 0; i < n; i++) {
    *p++ = ur[i];
    *p++ = ug[i];
    *p++ = ub[i];
  }

  delete[] ur;
  delete[] ug;
  delete[] ub;

  return rgb;
}

unsigned char *FitsDisplay::XYZ_AdobeRGB(size_t n, const float *X, 
					 const float *Y, const float *Z)
{
  // http://en.wikipedia.org/wiki/Adobe_RGB_color_space
  // http://www.adobe.com/digitalimag/pdfs/AdobeRGB1998.pdf

  float *r = new float[n];
  float *g = new float[n];
  float *b = new float[n];

  for(size_t i = 0; i < n; i++) {
    float x = X[i]/xtemp;
    float y = Y[i]/ytemp;
    float z = Z[i]/ztemp;

    r[i] =  2.04159f*x - 0.56501f*y - 0.34473f*z;
    g[i] = -0.96924f*x + 1.87597f*y + 0.04156f*z;
    b[i] =  0.01344f*x - 0.11836f*y + 1.01517f*z;
  }

  wxASSERT(cspace == "AdobeRGB");

  float *r1 = AdobeGamma(n,r);
  float *g1 = AdobeGamma(n,g);
  float *b1 = AdobeGamma(n,b);

  delete[] r;
  delete[] g;
  delete[] b;

  unsigned char *ur = 0;
  unsigned char *ug = 0;
  unsigned char *ub = 0;

  ur = Crange(n,r1);
  ug = Crange(n,g1);
  ub = Crange(n,b1);

  delete[] r1;
  delete[] g1;
  delete[] b1;


  unsigned char *rgb = (unsigned char *) malloc(3*n);
  unsigned char *p = rgb;

  for(size_t i = 0; i < n; i++) {
    *p++ = ur[i];
    *p++ = ug[i];
    *p++ = ub[i];
  }

  delete[] ur;
  delete[] ug;
  delete[] ub;

  return rgb;
}



// ----------  FitsHisto

class FitsHistoData : public wxObjectRefData
{
public:
  FitsHistoData();
  FitsHistoData(int, int *, float *);
  FitsHistoData(const FitsHistoData&);
  FitsHistoData& operator = (const FitsHistoData&);
  virtual ~FitsHistoData();

  int nbin;
  int *hist;
  float *cents;
};


FitsHistoData::FitsHistoData(): nbin(0),hist(0),cents(0) {}

FitsHistoData::FitsHistoData(int n, int *h, float *t): 
  nbin(n),hist(h),cents(t) {}


FitsHistoData::FitsHistoData(const FitsHistoData& copy) 
{ 
  wxFAIL_MSG("FitsHistoData ----- WE ARE REALY NEED COPY CONSTRUCTOR -----");
}

FitsHistoData& FitsHistoData::operator = (const FitsHistoData& other)
{
  wxFAIL_MSG("*** FitsHistoData: WE ARE REALLY NEED ASSIGNMENT CONSTRUCTOR");
  return *this;
}


FitsHistoData::~FitsHistoData() 
{
  delete[] hist;
  delete[] cents;
}


FitsHisto::FitsHisto(): wbin(0.0), xmin(0.0), xmax(0.0), npixels(0) {}

FitsHisto::FitsHisto(int n, int *h, float *c)
{
  UnRef();
  SetRefData(new FitsHistoData(n,h,c));
  xmin = c[0];
  xmax = c[n-1];
  wbin = (xmax - xmin)/2.0;
  npixels = 0;
}


FitsHisto::FitsHisto(const FitsArray& a, int nb)
{
  wxASSERT(a.IsOk());

  //nbin = int(log((double) array->Npixels())/log(2.0) + 1.0);
  // http://en.wikipedia.org/wiki/Histogram - Sturges' formula
  // * perhaps unusable due to low number of bins (21 for 1e6)

  int nbin = nb > 0 ? nb : 256;
  FitsArrayStat array(a);

  // bin width
  wbin = 3.5*array.Mad()/powf(float(array.Npixels()),0.333);

  xmin = array.Med() - nbin/2*wbin;
  xmax = array.Med() + nbin/2*wbin;
  npixels = array.Npixels();

  Create(nbin,array);
}

FitsHisto::FitsHisto(const FitsArray& a, float x0, float x1, int nb):
  xmin(x0), xmax(x1)
{
  wxASSERT(a.IsOk());
  int nbin = nb > 0 ? nb : 256;

  FitsArrayStat array(a);

  // bin width
  wbin = max((xmax - xmin)/float(nbin),float(10.0*FLT_EPSILON));
  npixels = array.Npixels();

  Create(nbin,array);
}

FitsHisto::~FitsHisto() {}

bool FitsHisto::IsOk() const
{
  const FitsHistoData *data = static_cast<const FitsHistoData *>(m_refData);
  return data && data->nbin > 0 && data->hist && data->cents;
}

void FitsHisto::Create(int nbin, const FitsArrayStat& a)
{
  wxASSERT(a.IsOk() && wbin > FLT_EPSILON);
  // http://en.wikipedia.org/wiki/Kernel_density_estimation ?

  int *hist = new int[nbin];
  float *cents = new float[nbin];

  for(int n = 0; n < nbin; n++) hist[n] = 0;

  for(int n = 0; n < nbin; n++) cents[n] = xmin + n*wbin;

  for(int i = 0; i < a.Npixels(); i=i+107) {
    int n = int((a.Pixel(i) - xmin)/wbin);
    if( 0 <= n  && n < nbin ) hist[n] += 107;
  }

  UnRef();
  SetRefData(new FitsHistoData(nbin,hist,cents));


#ifdef __WXDEBUG__
  wxFFileOutputStream file("/tmp/hist","w+");
  wxTextOutputStream cout(file);
  for(int i = 0; i < nbin; i++)
    cout << i << " " << cents[i] << " " << hist[i] << endl;
#endif
}

int FitsHisto::Hist(int n) const 
{ 
  const FitsHistoData *data = static_cast<const FitsHistoData *>(m_refData);
  wxASSERT(data && 0 <= n && n < data->nbin && data->hist);

  return data->hist[n];
}

float FitsHisto::Cents(int n) const 
{ 
  const FitsHistoData *data = static_cast<const FitsHistoData *>(m_refData);
  wxASSERT(data && 0 <= n && n < data->nbin && data->cents);
  return data->cents[n]; 
}

float FitsHisto::CentsMin() const 
{ 
  const FitsHistoData *data = static_cast<const FitsHistoData *>(m_refData);
  wxASSERT(data && data->nbin > 0 && data->cents && data->hist);

  for(int i = 0; i < data->nbin-1; i++)
    if( data->hist[i] > 0 )
      return data->cents[i];

  return data->cents[0];
}

float FitsHisto::CentsMax() const 
{ 
  const FitsHistoData *data = static_cast<const FitsHistoData *>(m_refData);
  wxASSERT(data && data->nbin > 0 && data->cents && data->hist);

  for(int i = data->nbin-1; i > 0; i--)
    if( data->hist[i] > 0 )
      return data->cents[i];

  return data->cents[data->nbin-1];
}

float FitsHisto::BinWidth() const 
{ 
  return wbin;
}


int FitsHisto::NBins() const 
{ 
  const FitsHistoData *data = static_cast<const FitsHistoData *>(m_refData);
  if( data )
    return data->nbin;
  else
    return 0;
}




// ---  FitsBitmap

// A simple, thread safe, replacement for wxImage.
// Memory should be managed via malloc/free functions.

class FitsBitmapData: public wxObjectRefData
{
public:
  FitsBitmapData();
  FitsBitmapData(int, int, unsigned char *);
  virtual ~FitsBitmapData();

  int width, height;
  unsigned char *rgb;
};
 

FitsBitmapData::FitsBitmapData(): width(0),height(0),rgb(0) {}

FitsBitmapData::FitsBitmapData(int w, int h, unsigned char *d): 
  width(w),height(h),rgb(d) {}


FitsBitmapData::~FitsBitmapData() 
{ 
  free(rgb);
}


FitsBitmap::FitsBitmap() {}

FitsBitmap::FitsBitmap(int w,int h, unsigned char *d)
{
  wxASSERT(w > 0 && h > 0 && d);
  UnRef();
  SetRefData(new FitsBitmapData(w,h,d));
}

FitsBitmap::~FitsBitmap() {}

int FitsBitmap::GetWidth() const
{
  FitsBitmapData *data = static_cast<FitsBitmapData *>(m_refData);
  wxASSERT(data);
  return data->width;
}

int FitsBitmap::GetHeight() const
{
  FitsBitmapData *data = static_cast<FitsBitmapData *>(m_refData);
  wxASSERT(data);
  return data->height;
}

const unsigned char *FitsBitmap::GetRGB() const
{
  FitsBitmapData *data = static_cast<FitsBitmapData *>(m_refData);
  wxASSERT(data);
  return data->rgb;
}

unsigned char *FitsBitmap::NewRGB() const
{
  FitsBitmapData *data = static_cast<FitsBitmapData *>(m_refData);
  wxASSERT(data);
  int npix = 3*data->width*data->height;
  unsigned char *rgb = (unsigned char *) malloc(npix);
  memcpy(rgb,data->rgb,npix);
  return rgb;
}

unsigned char *FitsBitmap::NewTopsyTurvyRGB() const
{
  FitsBitmapData *data = static_cast<FitsBitmapData *>(m_refData);
  wxASSERT(data);
  int npix = 3*data->width;
  unsigned char *rgb = (unsigned char *) malloc(npix*data->height);
  for(int i = 0; i < data->height; i++)
    memcpy(rgb+i*npix,data->rgb+(data->height - i - 1)*npix,npix);
  return rgb;
}

bool FitsBitmap::IsOk() const
{
  FitsBitmapData *data = static_cast<FitsBitmapData *>(m_refData);
  return data && data->width > 0 && data->height > 0 && data->rgb;
}
